library(ggplot2)
ll.nbinom <- function(par, x, d) {
if (min(par) < 0) {
return(Inf)
}
-sum(dnbinom(x, par[1], (par[2]/d)/((par[2]/d) + 1), log=TRUE))
}
lk.nbinom <- function(par, x, d, priors=list(a=c(1, 0.5), b=c(1, 0.5))) {
if (min(par) < 0) {
return(-Inf)
}
r <- par[1]
p <- par[2]/(par[2] + d)
sum(dnbinom(x, r, p, log=TRUE)) + dgamma(par[1], priors$a[1], priors$a[2], log=TRUE) + dgamma(par[2], priors$b[1], priors$b[2], log=TRUE)
}
rtnorm <- function(n, mean = 0, sd = 1, lower=-Inf, upper=Inf) {
# range is a vector of two values
F.a <- pnorm(lower, mean = mean, sd = sd)
F.b <- pnorm(upper, mean = mean, sd = sd)
u <- runif(n, min = F.a, max = F.b)
qnorm(u, mean = mean, sd = sd)
}
metropolis <- function(target, initial, iterations=10000, warmup=2000, thin=100) {
k <- length(initial)
chain <- matrix(0, iterations + warmup, k)
accepted <- 0
chain[1, ] <- initial
for (i in 2:(iterations + warmup)) {
current <- chain[i - 1, ]
candidate <- rnorm(k, mean=current, sd=1)
# candidate <- rtnorm(k, mean=current, sd=1, lower=0, upper=)
A <- exp(target(candidate) - target(current))
if (runif(1) < A) {
accepted <- accepted + 1
chain[i, ] <- candidate
} else {
chain[i, ] <- current
}
}
# print(accepted/(iterations+warmup))
return(chain[seq(warmup, iterations + warmup, thin), ])
}
gibbs <- function(target, initial, iterations=10000, warmup=2000, thin=100) {
k <- length(initial)
chain <- matrix(0, iterations + warmup, k)
accepted <- 0
chain[1, ] <- initial
for (i in 2:(iterations + warmup)) {
current <- chain[i - 1, ]
candidate <- chain[i - 1, ]
for (j in 1:k) {
candidate[j] <- rnorm(1, mean=current[j], sd=1)
# candidate[j] <- rtnorm(1, mean=current[j], sd=1, lower=0, upper=Inf)
A <- exp(target(candidate) - target(current))
if (runif(1) < A) {
current <- candidate
} else {
candidate <- current
}
}
chain[i, ] <- current
}
# print(accepted/(iterations+warmup))
return(chain[seq(warmup, iterations + warmup, thin), ])
}
plot(density(rgamma(1000, 1, 1)))
lines(density(rgamma(1000, 1, 0.1)), col=2)

n <- 100
iterations <- 100
cors_true <- rep(0, iterations)
cors_est <- rep(0, iterations)
cors_est_lb <- rep(0, iterations)
cors_est_ub <- rep(0, iterations)
mean_sampling <- rep(0, iterations)
bad_chains <- rep(0, iterations)
i <- 1
for (i in 1:iterations) {
print(i)
a <- runif(1, 1, 10)
b <- runif(1, 1, 10)
# Do variable sampling effort.
mean_sampling[i] <- runif(1, 1, 20)
d <- rpois(n, mean_sampling[i])
d[d == 0] = 1
alpha <- rgamma(n, a, b)
alpha_hat <- rpois(n, alpha*d)/d
cors_true[i] <- cor(alpha, alpha_hat)
# cors_true[i] <- var(alpha)/var(alpha_hat)
# parameters <- optim(c(1, 1), ll.nbinom, x=as.integer(alpha_hat * d), d=d)$par
# a_ <- parameters[1]
# b_ <- parameters[2]
target <- function(par) lk.nbinom(par, as.integer(alpha_hat * d), d)
chain <- gibbs(target, c(1, 1), warmup=2000, iterations=20000, thin=10)
plot(chain[, 1], type="l", col=1, ylim=c(0, max(chain)))
lines(chain[, 2], type="l", col=2)
qt <- quantile(chain[, 2], probs=c(0.025, 0.5, 0.975))
b_ <- qt[2]
b_lb <- qt[1]
b_ub <- qt[3]
cors_est[i] <- sqrt(n/(n + b_ * sum(d^-1)))
cors_est_ub[i] <- sqrt(n/(n + b_ub * sum(d^-1)))
cors_est_lb[i] <- sqrt(n/(n + b_lb * sum(d^-1)))
}
[1] 1
[1] 2

[1] 3

[1] 4

[1] 5

[1] 6

[1] 7

[1] 8

[1] 9

[1] 10

[1] 11

[1] 12

[1] 13

[1] 14

[1] 15

[1] 16

[1] 17

[1] 18

[1] 19

[1] 20

[1] 21

[1] 22

[1] 23

[1] 24

[1] 25

[1] 26

[1] 27

[1] 28

[1] 29

[1] 30

[1] 31

[1] 32

[1] 33

[1] 34

[1] 35

[1] 36

[1] 37

[1] 38

[1] 39

[1] 40

[1] 41

[1] 42

[1] 43

[1] 44

[1] 45

[1] 46

[1] 47

[1] 48

[1] 49

[1] 50

[1] 51

[1] 52

[1] 53

[1] 54

[1] 55

[1] 56

[1] 57

[1] 58

[1] 59

[1] 60

[1] 61

[1] 62

[1] 63

[1] 64

[1] 65

[1] 66

[1] 67

[1] 68

[1] 69

[1] 70

[1] 71

[1] 72

[1] 73

[1] 74

[1] 75

[1] 76

[1] 77

[1] 78

[1] 79

[1] 80

[1] 81

[1] 82

[1] 83

[1] 84

[1] 85

[1] 86

[1] 87

[1] 88

[1] 89

[1] 90

[1] 91

[1] 92

[1] 93

[1] 94

[1] 95

[1] 96

[1] 97

[1] 98

[1] 99

[1] 100


qplot(cors_true, cors_est) +
geom_errorbar(aes(ymin=cors_est_lb, ymax=cors_est_ub)) +
geom_smooth() +
geom_abline() +
xlab("True correlation") +
ylab("Estimated correlation")

qplot(mean_sampling, abs(cors_true - cors_est)) +
# geom_errorbar(aes(ymin=abs(cors_true-cors_est_lb), ymax=abs(cors_true-cors_est_ub))) +
geom_smooth()

a <- runif(1, 1, 10)
b <- runif(1, 1, 10)
X <- seq(0.1, 100, 0.1)
alpha <- rgamma(n, a, b)
priors <- list(a=c(1, 1), b=c(1, 1))
d <- rep(1, n)
alpha_hat <- rpois(n, alpha*d)/d
plot(X, sapply(X, function(x) lk.nbinom(c(a, x), as.integer(alpha_hat * d), d, priors=priors)), type="l", col=1)
d <- rep(10, n)
alpha_hat <- rpois(n, alpha*d)/d
lines(X, sapply(X, function(x) lk.nbinom(c(a, x), as.integer(alpha_hat * d), d, priors=priors)), type="l", col=2)
d <- rep(100, n)
alpha_hat <- rpois(n, alpha*d)/d
lines(X, sapply(X, function(x) lk.nbinom(c(a, x), as.integer(alpha_hat * d), d, priors=priors)), type="l", col=3)


LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CmxpYnJhcnkoZ2dwbG90MikKYGBgCgpgYGB7cn0KbGwubmJpbm9tIDwtIGZ1bmN0aW9uKHBhciwgeCwgZCkgewogIGlmIChtaW4ocGFyKSA8IDApIHsKICAgIHJldHVybihJbmYpCiAgfQogIC1zdW0oZG5iaW5vbSh4LCBwYXJbMV0sIChwYXJbMl0vZCkvKChwYXJbMl0vZCkgKyAxKSwgbG9nPVRSVUUpKQp9Cgpsay5uYmlub20gPC0gZnVuY3Rpb24ocGFyLCB4LCBkLCBwcmlvcnM9bGlzdChhPWMoMSwgMC41KSwgYj1jKDEsIDAuNSkpKSB7CiAgaWYgKG1pbihwYXIpIDwgMCkgewogICAgcmV0dXJuKC1JbmYpCiAgfQogIHIgPC0gcGFyWzFdCiAgcCA8LSBwYXJbMl0vKHBhclsyXSArIGQpCiAgc3VtKGRuYmlub20oeCwgciwgcCwgbG9nPVRSVUUpKSArIGRnYW1tYShwYXJbMV0sIHByaW9ycyRhWzFdLCBwcmlvcnMkYVsyXSwgbG9nPVRSVUUpICsgZGdhbW1hKHBhclsyXSwgcHJpb3JzJGJbMV0sIHByaW9ycyRiWzJdLCBsb2c9VFJVRSkKfQoKcnRub3JtIDwtIGZ1bmN0aW9uKG4sIG1lYW4gPSAwLCBzZCA9IDEsIGxvd2VyPS1JbmYsIHVwcGVyPUluZikgewogICMgcmFuZ2UgaXMgYSB2ZWN0b3Igb2YgdHdvIHZhbHVlcwogIAogIEYuYSA8LSBwbm9ybShsb3dlciwgbWVhbiA9IG1lYW4sIHNkID0gc2QpCiAgRi5iIDwtIHBub3JtKHVwcGVyLCBtZWFuID0gbWVhbiwgc2QgPSBzZCkKICAKICB1IDwtIHJ1bmlmKG4sIG1pbiA9IEYuYSwgbWF4ID0gRi5iKQogIAogIHFub3JtKHUsIG1lYW4gPSBtZWFuLCBzZCA9IHNkKQp9CgptZXRyb3BvbGlzIDwtIGZ1bmN0aW9uKHRhcmdldCwgaW5pdGlhbCwgaXRlcmF0aW9ucz0xMDAwMCwgd2FybXVwPTIwMDAsIHRoaW49MTAwKSB7CiAgayA8LSBsZW5ndGgoaW5pdGlhbCkKICBjaGFpbiA8LSBtYXRyaXgoMCwgaXRlcmF0aW9ucyArIHdhcm11cCwgaykKICBhY2NlcHRlZCA8LSAwCiAgCiAgY2hhaW5bMSwgXSA8LSBpbml0aWFsCiAgCiAgZm9yIChpIGluIDI6KGl0ZXJhdGlvbnMgKyB3YXJtdXApKSB7CiAgICBjdXJyZW50IDwtIGNoYWluW2kgLSAxLCBdCiAgICBjYW5kaWRhdGUgPC0gcm5vcm0oaywgbWVhbj1jdXJyZW50LCBzZD0xKQogICAgIyBjYW5kaWRhdGUgPC0gcnRub3JtKGssIG1lYW49Y3VycmVudCwgc2Q9MSwgbG93ZXI9MCwgdXBwZXI9KQogICAgQSA8LSBleHAodGFyZ2V0KGNhbmRpZGF0ZSkgLSB0YXJnZXQoY3VycmVudCkpCiAgICBpZiAocnVuaWYoMSkgPCBBKSB7CiAgICAgIGFjY2VwdGVkIDwtIGFjY2VwdGVkICsgMQogICAgICBjaGFpbltpLCBdIDwtIGNhbmRpZGF0ZQogICAgfSBlbHNlIHsKICAgICAgY2hhaW5baSwgXSA8LSBjdXJyZW50CiAgICB9CiAgfQogICMgcHJpbnQoYWNjZXB0ZWQvKGl0ZXJhdGlvbnMrd2FybXVwKSkKICByZXR1cm4oY2hhaW5bc2VxKHdhcm11cCwgaXRlcmF0aW9ucyArIHdhcm11cCwgdGhpbiksIF0pCn0KCmdpYmJzIDwtIGZ1bmN0aW9uKHRhcmdldCwgaW5pdGlhbCwgaXRlcmF0aW9ucz0xMDAwMCwgd2FybXVwPTIwMDAsIHRoaW49MTAwKSB7CiAgayA8LSBsZW5ndGgoaW5pdGlhbCkKICBjaGFpbiA8LSBtYXRyaXgoMCwgaXRlcmF0aW9ucyArIHdhcm11cCwgaykKICBhY2NlcHRlZCA8LSAwCiAgCiAgY2hhaW5bMSwgXSA8LSBpbml0aWFsCiAgCiAgZm9yIChpIGluIDI6KGl0ZXJhdGlvbnMgKyB3YXJtdXApKSB7CiAgICBjdXJyZW50IDwtIGNoYWluW2kgLSAxLCBdCiAgICBjYW5kaWRhdGUgPC0gY2hhaW5baSAtIDEsIF0KICAgIGZvciAoaiBpbiAxOmspIHsKICAgICAgY2FuZGlkYXRlW2pdIDwtIHJub3JtKDEsIG1lYW49Y3VycmVudFtqXSwgc2Q9MSkKICAgICAgIyBjYW5kaWRhdGVbal0gPC0gcnRub3JtKDEsIG1lYW49Y3VycmVudFtqXSwgc2Q9MSwgbG93ZXI9MCwgdXBwZXI9SW5mKQogICAgICBBIDwtIGV4cCh0YXJnZXQoY2FuZGlkYXRlKSAtIHRhcmdldChjdXJyZW50KSkKICAgICAgaWYgKHJ1bmlmKDEpIDwgQSkgewogICAgICAgIGN1cnJlbnQgPC0gY2FuZGlkYXRlCiAgICAgIH0gZWxzZSB7CiAgICAgICAgY2FuZGlkYXRlIDwtIGN1cnJlbnQKICAgICAgfQogICAgfQogICAgY2hhaW5baSwgXSA8LSBjdXJyZW50CiAgfQogICMgcHJpbnQoYWNjZXB0ZWQvKGl0ZXJhdGlvbnMrd2FybXVwKSkKICByZXR1cm4oY2hhaW5bc2VxKHdhcm11cCwgaXRlcmF0aW9ucyArIHdhcm11cCwgdGhpbiksIF0pCn0KYGBgCgpgYGB7cn0KcGxvdChkZW5zaXR5KHJnYW1tYSgxMDAwLCAxLCAxKSkpCmxpbmVzKGRlbnNpdHkocmdhbW1hKDEwMDAsIDEsIDAuMSkpLCBjb2w9MikKYGBgCgpgYGB7cn0KbiA8LSAxMDAKaXRlcmF0aW9ucyA8LSAxMDAKY29yc190cnVlIDwtIHJlcCgwLCBpdGVyYXRpb25zKQpjb3JzX2VzdCA8LSByZXAoMCwgaXRlcmF0aW9ucykKY29yc19lc3RfbGIgPC0gcmVwKDAsIGl0ZXJhdGlvbnMpCmNvcnNfZXN0X3ViIDwtIHJlcCgwLCBpdGVyYXRpb25zKQptZWFuX3NhbXBsaW5nIDwtIHJlcCgwLCBpdGVyYXRpb25zKQpiYWRfY2hhaW5zIDwtIHJlcCgwLCBpdGVyYXRpb25zKQoKaSA8LSAxCgpmb3IgKGkgaW4gMTppdGVyYXRpb25zKSB7CiAgcHJpbnQoaSkKICBhIDwtIHJ1bmlmKDEsIDEsIDEwKQogIGIgPC0gcnVuaWYoMSwgMSwgMTApCiAgCiAgIyBEbyB2YXJpYWJsZSBzYW1wbGluZyBlZmZvcnQuCiAgbWVhbl9zYW1wbGluZ1tpXSA8LSBydW5pZigxLCAxLCAyMCkKICBkIDwtIHJwb2lzKG4sIG1lYW5fc2FtcGxpbmdbaV0pCiAgZFtkID09IDBdID0gMQogIAogIGFscGhhIDwtIHJnYW1tYShuLCBhLCBiKQogIGFscGhhX2hhdCA8LSBycG9pcyhuLCBhbHBoYSpkKS9kCiAgCiAgY29yc190cnVlW2ldIDwtIGNvcihhbHBoYSwgYWxwaGFfaGF0KQogICMgY29yc190cnVlW2ldIDwtIHZhcihhbHBoYSkvdmFyKGFscGhhX2hhdCkKCiAgIyBwYXJhbWV0ZXJzIDwtIG9wdGltKGMoMSwgMSksIGxsLm5iaW5vbSwgeD1hcy5pbnRlZ2VyKGFscGhhX2hhdCAqIGQpLCBkPWQpJHBhcgogICMgYV8gPC0gcGFyYW1ldGVyc1sxXQogICMgYl8gPC0gcGFyYW1ldGVyc1syXQogIAogIHRhcmdldCA8LSBmdW5jdGlvbihwYXIpIGxrLm5iaW5vbShwYXIsIGFzLmludGVnZXIoYWxwaGFfaGF0ICogZCksIGQpCiAgY2hhaW4gPC0gZ2liYnModGFyZ2V0LCBjKDEsIDEpLCB3YXJtdXA9MjAwMCwgaXRlcmF0aW9ucz0yMDAwMCwgdGhpbj0xMCkKCiAgcGxvdChjaGFpblssIDFdLCB0eXBlPSJsIiwgY29sPTEsIHlsaW09YygwLCBtYXgoY2hhaW4pKSkKICBsaW5lcyhjaGFpblssIDJdLCB0eXBlPSJsIiwgY29sPTIpCiAgICAKICBxdCA8LSBxdWFudGlsZShjaGFpblssIDJdLCBwcm9icz1jKDAuMDI1LCAwLjUsIDAuOTc1KSkKCiAgYl8gPC0gcXRbMl0KICBiX2xiIDwtIHF0WzFdCiAgYl91YiA8LSBxdFszXQogIAogIGNvcnNfZXN0W2ldIDwtIHNxcnQobi8obiArIGJfICogc3VtKGReLTEpKSkKICBjb3JzX2VzdF91YltpXSA8LSBzcXJ0KG4vKG4gKyBiX3ViICogc3VtKGReLTEpKSkKICBjb3JzX2VzdF9sYltpXSA8LSBzcXJ0KG4vKG4gKyBiX2xiICogc3VtKGReLTEpKSkKfQoKcXBsb3QoY29yc190cnVlLCBjb3JzX2VzdCkgKwogIGdlb21fZXJyb3JiYXIoYWVzKHltaW49Y29yc19lc3RfbGIsIHltYXg9Y29yc19lc3RfdWIpKSArCiAgZ2VvbV9zbW9vdGgoKSArCiAgZ2VvbV9hYmxpbmUoKSArCiAgeGxhYigiVHJ1ZSBjb3JyZWxhdGlvbiIpICsKICB5bGFiKCJFc3RpbWF0ZWQgY29ycmVsYXRpb24iKQoKcXBsb3QobWVhbl9zYW1wbGluZywgYWJzKGNvcnNfdHJ1ZSAtIGNvcnNfZXN0KSkgKwogICMgZ2VvbV9lcnJvcmJhcihhZXMoeW1pbj1hYnMoY29yc190cnVlLWNvcnNfZXN0X2xiKSwgeW1heD1hYnMoY29yc190cnVlLWNvcnNfZXN0X3ViKSkpICsKICBnZW9tX3Ntb290aCgpCmBgYAoKYGBge3J9CmEgPC0gcnVuaWYoMSwgMSwgMTApCmIgPC0gcnVuaWYoMSwgMSwgMTApCgpYIDwtIHNlcSgwLjEsIDEwMCwgMC4xKQoKYWxwaGEgPC0gcmdhbW1hKG4sIGEsIGIpCgpwcmlvcnMgPC0gbGlzdChhPWMoMSwgMSksIGI9YygxLCAxKSkKCmQgPC0gcmVwKDEsIG4pCmFscGhhX2hhdCA8LSBycG9pcyhuLCBhbHBoYSpkKS9kCnBsb3QoWCwgc2FwcGx5KFgsIGZ1bmN0aW9uKHgpIGxrLm5iaW5vbShjKGEsIHgpLCBhcy5pbnRlZ2VyKGFscGhhX2hhdCAqIGQpLCBkLCBwcmlvcnM9cHJpb3JzKSksIHR5cGU9ImwiLCBjb2w9MSkKCmQgPC0gcmVwKDEwLCBuKQphbHBoYV9oYXQgPC0gcnBvaXMobiwgYWxwaGEqZCkvZApsaW5lcyhYLCBzYXBwbHkoWCwgZnVuY3Rpb24oeCkgbGsubmJpbm9tKGMoYSwgeCksIGFzLmludGVnZXIoYWxwaGFfaGF0ICogZCksIGQsIHByaW9ycz1wcmlvcnMpKSwgdHlwZT0ibCIsIGNvbD0yKQoKZCA8LSByZXAoMTAwLCBuKQphbHBoYV9oYXQgPC0gcnBvaXMobiwgYWxwaGEqZCkvZApsaW5lcyhYLCBzYXBwbHkoWCwgZnVuY3Rpb24oeCkgbGsubmJpbm9tKGMoYSwgeCksIGFzLmludGVnZXIoYWxwaGFfaGF0ICogZCksIGQsIHByaW9ycz1wcmlvcnMpKSwgdHlwZT0ibCIsIGNvbD0zKQpgYGAKCmBgYHtyfQpuIDwtIDEwMAppdGVyYXRpb25zIDwtIDEwMAoKbWVhbl9zYW1wbGluZyA8LSByZXAoMCwgaXRlcmF0aW9ucykKdGVzdF9zdGF0aXN0aWNzIDwtIHJlcCgwLCBpdGVyYXRpb25zKQoKaSA8LSAxCgpmb3IgKGkgaW4gMTppdGVyYXRpb25zKSB7CiAgYSA8LSBydW5pZigxLCAxLCAxMCkKICBiIDwtIHJ1bmlmKDEsIDEsIDEwKQogIAogICMgRG8gdmFyaWFibGUgc2FtcGxpbmcgZWZmb3J0LgogIG1lYW5fc2FtcGxpbmdbaV0gPC0gcnVuaWYoMSwgMSwgMjApCiAgZCA8LSBycG9pcyhuLCBtZWFuX3NhbXBsaW5nW2ldKQogIGRbZCA9PSAwXSA9IDEKICAKICBhbHBoYSA8LSByZ2FtbWEobiwgYSwgYikKICBhbHBoYV9oYXQgPC0gcnBvaXMobiwgYWxwaGEqZCkvZAogIGluZGV4X2RpZmYgPC0gYWxwaGEgLSBhbHBoYV9oYXQKICAjIGNvcihpbmRleF9kaWZmLCBybm9ybShsZW5ndGgoaW5kZXhfZGlmZiksIG1lYW4oaW5kZXhfZGlmZiksIHNkKGluZGV4X2RpZmYpKSkKICBxcS55IDwtIChhcy52ZWN0b3IocXVhbnRpbGUoaW5kZXhfZGlmZiwgcHJvYnM9c2VxKDAsIDEsIDAuMDEpKSkgLSBtZWFuKGluZGV4X2RpZmYpKS9zZChpbmRleF9kaWZmKQogIHFxLnggPC0gYXMudmVjdG9yKHF1YW50aWxlKHJub3JtKDFlNCwgMCwgMSksIHByb2JzPXNlcSgwLCAxLCAwLjAxKSkpCiAgCiAgIyBxcGxvdChxcS54LCBxcS55KSArCiAgIyAgIGdlb21fYWJsaW5lKCkgKwogICMgICBnZW9tX3Ntb290aChtZXRob2Q9ImxtIikKICAKICB0ZXN0X3N0YXRpc3RpY3NbaV0gPC0gc3VtKGFicyhxcS54IC0gcXEueSkpCgogICMgc3cgPC0gc2hhcGlyby50ZXN0KGFscGhhIC0gYWxwaGFfaGF0KQogICMgdGVzdF9zdGF0aXN0aWNzW2ldIDwtIHN3JHN0YXRpc3RpY1tbMV1dCn0KCnFwbG90KG1lYW5fc2FtcGxpbmcsIHRlc3Rfc3RhdGlzdGljcykgKwogICMgZ2VvbV9lcnJvcmJhcihhZXMoeW1pbj1hYnMoY29yc190cnVlLWNvcnNfZXN0X2xiKSwgeW1heD1hYnMoY29yc190cnVlLWNvcnNfZXN0X3ViKSkpICsKICBnZW9tX3Ntb290aCgpICsKICB4bGFiKCJNZWFuIHNhbXBsaW5nIChuKSIpICsKICB5bGFiKCJRUSBub3JtYWxpdHkgZXJyb3IiKQpgYGAK